/*


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-lg93.c:

   LG93 functional: gradient correction GGA exchange functional, 
   with the full LG93 exchange actually SLATER + LG93x.

   Implemented by David Wilson (david.wilson@latrobe.edu.au) Jun 2005

   Automatically generated code implementing LG93 functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional lg93Functional;" to 'functionals.h'
    2. add "&lg93Functional," to 'functionals.c'
    3. add "fun-lg93.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

ad:   1.0*10^(-8);
a4:   29.790;
a6:   22.417;
a8:   12.119;
a10:  1570.1;
a12:  55.944;
b:    0.024974;
C0:   (ad+0.1234)/b;

xa:   grada/(rhoa^(4/3));
xb:   gradb/(rhob^(4/3));

sa:   0.5*((3*%PI^2)^(-1/3))*xa;
sb:   0.5*((3*%PI^2)^(-1/3))*xb;

F1a:    (1+C0*sa^2+a4*sa^4+a6*sa^6+a8*sa^8+a10*sa^10+a12*sa^12)^b;
F1b:    (1+C0*sb^2+a4*sb^4+a6*sb^6+a8*sb^8+a10*sb^10+a12*sb^12)^b;
F2a:    1+ad*sa^2;
F2b:    1+ad*sb^2;
Exa:    rhoa^(4/3)*F1a/F2a;
Exb:    rhob^(4/3)*F1b/F2b;

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer lg93_isgga(void) { return 1; } /* FIXME: detect! */
static integer lg93_read(const char *conf_line);
static real lg93_energy(const FunDensProp* dp);
static void lg93_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void lg93_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void lg93_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void lg93_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional LG93xFunctional = {
  "LG93x",       /* name */
  lg93_isgga,   /* gga-corrected */
   1,
  lg93_read,
  NULL,
  lg93_energy,
  lg93_first,
  lg93_second,
  lg93_third,
  lg93_fourth
};
 
/* IMPLEMENTATION PART */
static integer
lg93_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
lg93_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12;

    t1 = 1/pow(3.0,0.666666666666667);
    t2 = 1/pow(3.141592653589793,1.333333333333333);
    t3 = pow(grada,2.0);
    t4 = 1/pow(rhoa,2.666666666666667);
    t5 = 1/pow(3.141592653589793,8.0);
    t6 = 1/pow(3.0,0.333333333333333);
    t7 = 1/pow(3.141592653589793,6.666666666666667);
    t8 = 1/pow(3.141592653589793,5.333333333333333);
    t9 = 1/pow(3.141592653589793,4.0);
    t10 = 1/pow(3.141592653589793,2.666666666666667);
    t11 = pow(gradb,2.0);
    t12 = 1/pow(rhob,2.666666666666667);

   /* code */
    res = pow(rhob,1.333333333333333)*pow(1.6861979166666666E-4*
        t5*pow(gradb,12.0)/pow(rhob,16.0)+0.056788917824074*t6*t7*
        pow(gradb,10.0)/pow(rhob,13.33333333333333)+0.005259982638889*
        t1*t8*pow(gradb,8.0)/pow(rhob,10.66666666666667)+0.038918402777778*
        t9*pow(gradb,6.0)/pow(rhob,8.0)+0.620625*t10*t6*pow(gradb,
        4.0)/pow(rhob,5.333333333333333)+1.235284796188036*t1*t2*t11*
        t12+1.0,0.024974)/(2.5E-9*t1*t2*t11*t12+1.0)+pow(rhoa,1.333333333333333)*
        pow(1.6861979166666666E-4*t5*pow(grada,12.0)/pow(rhoa,16.0)+
        0.056788917824074*t6*t7*pow(grada,10.0)/pow(rhoa,13.33333333333333)+
        0.005259982638889*t1*t8*pow(grada,8.0)/pow(rhoa,10.66666666666667)+
        0.038918402777778*t9*pow(grada,6.0)/pow(rhoa,8.0)+0.620625*
        t10*t6*pow(grada,4.0)/pow(rhoa,5.333333333333333)+1.235284796188036*
        t1*t2*t3*t4+1.0,0.024974)/(2.5E-9*t1*t2*t3*t4+1.0);

    return res;
}

static void
lg93_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27;

    t1 = 1/pow(3.0,0.666666666666667);
    t2 = 1/pow(3.141592653589793,1.333333333333333);
    t3 = pow(grada,2.0);
    t4 = 1/pow(rhoa,2.666666666666667);
    t5 = 2.5E-9*t1*t2*t3*t4+1.0;
    t6 = 1/pow(t5,2.0);
    t7 = 1/pow(3.141592653589793,8.0);
    t8 = pow(grada,12.0);
    t9 = 1/pow(rhoa,16.0);
    t10 = 1/pow(3.0,0.333333333333333);
    t11 = 1/pow(3.141592653589793,6.666666666666667);
    t12 = pow(grada,10.0);
    t13 = 1/pow(rhoa,13.33333333333333);
    t14 = 1/pow(3.141592653589793,5.333333333333333);
    t15 = pow(grada,8.0);
    t16 = 1/pow(rhoa,10.66666666666667);
    t17 = 1/pow(3.141592653589793,4.0);
    t18 = pow(grada,6.0);
    t19 = 1/pow(rhoa,8.0);
    t20 = 1/pow(3.141592653589793,2.666666666666667);
    t21 = pow(grada,4.0);
    t22 = 1/pow(rhoa,5.333333333333333);
    t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
        t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
        t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
        t10*t11*t12*t13+1.0;
    t24 = pow(t23,0.024974);
    t25 = 1/t5;
    t26 = 1/pow(t23,0.975026);
    t27 = pow(rhoa,1.333333333333333);

   /* code */
    res[0] = 0.024974*t25*t26*t27*(-0.002697916666667*t7*
        t8/pow(rhoa,17.0)-0.757185570987654*t10*t11*t12/pow(rhoa,14.33333333333333)-
        0.056106481481481*t1*t14*t15/pow(rhoa,11.66666666666667)-0.311347222222222*
        t17*t18/pow(rhoa,9.0)-3.31*t10*t20*t21/pow(rhoa,6.333333333333333)-
        3.294092789834761*t1*t2*t3/pow(rhoa,3.666666666666667))+6.666666666666667E-9*
        t1*t2*t24*t3*t6/pow(rhoa,2.333333333333333)+1.333333333333333*
        t24*t25*pow(rhoa,0.333333333333333);
    res[1] = 0.024974*t25*t26*t27*(0.0020234375*t7*t9*pow(grada,
        11.0)+0.567889178240741*t10*t11*t13*pow(grada,9.0)+0.042079861111111*
        t1*t14*t16*pow(grada,7.0)+0.233510416666667*t17*t19*pow(grada,
        5.0)+2.4825*t10*t20*t22*pow(grada,3.0)+2.470569592376071*t1*
        t2*grada*t4)-5.0E-9*t1*t2*t24*t6*grada/t27;
}

static void
lg93_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    lg93_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        lg93_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
lg93_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45;

    t1 = 1/pow(3.0,0.666666666666667);
    t2 = 1/pow(3.141592653589793,1.333333333333333);
    t3 = pow(grada,2.0);
    t4 = 1/pow(rhoa,2.666666666666667);
    t5 = 2.5E-9*t1*t2*t3*t4+1.0;
    t6 = 1/pow(t5,2.0);
    t7 = 1/pow(3.141592653589793,8.0);
    t8 = pow(grada,12.0);
    t9 = 1/pow(rhoa,16.0);
    t10 = 1/pow(3.0,0.333333333333333);
    t11 = 1/pow(3.141592653589793,6.666666666666667);
    t12 = pow(grada,10.0);
    t13 = 1/pow(rhoa,13.33333333333333);
    t14 = 1/pow(3.141592653589793,5.333333333333333);
    t15 = pow(grada,8.0);
    t16 = 1/pow(rhoa,10.66666666666667);
    t17 = 1/pow(3.141592653589793,4.0);
    t18 = pow(grada,6.0);
    t19 = 1/pow(rhoa,8.0);
    t20 = 1/pow(3.141592653589793,2.666666666666667);
    t21 = pow(grada,4.0);
    t22 = 1/pow(rhoa,5.333333333333333);
    t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
        t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
        t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
        t10*t11*t12*t13+1.0;
    t24 = pow(t23,0.024974);
    t25 = 1/pow(rhoa,2.333333333333333);
    t26 = 1/t5;
    t27 = pow(rhoa,0.333333333333333);
    t28 = 1/pow(rhoa,17.0);
    t29 = 1/pow(rhoa,14.33333333333333);
    t30 = 1/pow(rhoa,11.66666666666667);
    t31 = 1/pow(rhoa,9.0);
    t32 = 1/pow(rhoa,6.333333333333333);
    t33 = 1/pow(rhoa,3.666666666666667);
    t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
        t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
        t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
        t7*t8*t28;
    t35 = 1/pow(t23,0.975026);
    t36 = pow(rhoa,1.333333333333333);
    t37 = 1/t36;
    t38 = pow(grada,11.0);
    t39 = pow(grada,9.0);
    t40 = pow(grada,7.0);
    t41 = pow(grada,5.0);
    t42 = pow(grada,3.0);
    t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
        t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
        t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
        t7*t38*t9;
    t44 = 1/pow(t5,3.0);
    t45 = 1/pow(t23,1.975026);

   /* code */
    res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
        t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
    res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
        t6*t24*t37;
    res[2] = 0.024974*t26*t35*t36*(0.045864583333333*t7*t8/
        pow(rhoa,18.0)+10.85299318415638*t10*t11*t12/pow(rhoa,15.33333333333333)+
        0.65457561728395*t1*t14*t15/pow(rhoa,12.66666666666667)+2.802125*
        t17*t18/pow(rhoa,10.0)+20.96333333333333*t10*t20*t21/pow(rhoa,
        7.333333333333333)+12.07834022939412*t1*t2*t3/pow(rhoa,4.666666666666667))+
        2.96296296296296E-17*t10*t20*t21*t24*t44/pow(rhoa,6.0)-6.666666666666668E-9*
        t1*t2*t24*t3*t6/pow(rhoa,3.333333333333333)+0.444444444444444*
        t24*t26/pow(rhoa,0.666666666666667)-0.024350299324*t26*pow(t34,
        2.0)*t36*t45+0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*
        t1*t2*t3*t34*t6*t35*t25;
    res[3] = -2.222222222222222E-17*t10*t20*t24*t42*t44/pow(rhoa,
        5.0)-1.2487E-10*t1*t2*grada*t34*t6*t35*t37-0.024350299324*
        t34*t43*t26*t45*t36+0.024974*(-6.588185579669522*t1*t2*grada*
        t33-13.24*t10*t20*t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*
        t1*t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*
        t7*t38*t28)*t26*t35*t36+0.033298666666667*t43*t26*t35*t27+
        1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25+6.666666666666667E-9*
        t1*t2*grada*t6*t24*t25;
    res[4] = 1.666666666666667E-17*t10*t20*t24*t3*t44/pow(rhoa,
        4.0)-0.024350299324*t26*t36*pow(t43,2.0)*t45-2.4974E-10*t1*
        t2*grada*t43*t6*t35*t37-5.0E-9*t1*t2*t6*t24*t37+0.024974*(2.470569592376071*
        t1*t2*t4+7.4475*t10*t20*t3*t22+1.167552083333333*t17*t21*t19+
        0.294559027777778*t1*t14*t18*t16+5.111002604166666*t10*t11*
        t15*t13+0.0222578125*t7*t12*t9)*t26*t35*t36;

}

static void
lg93_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    lg93_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        lg93_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
lg93_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63;

    t1 = 1/pow(3.0,0.666666666666667);
    t2 = 1/pow(3.141592653589793,1.333333333333333);
    t3 = pow(grada,2.0);
    t4 = 1/pow(rhoa,2.666666666666667);
    t5 = 2.5E-9*t1*t2*t3*t4+1.0;
    t6 = 1/pow(t5,2.0);
    t7 = 1/pow(3.141592653589793,8.0);
    t8 = pow(grada,12.0);
    t9 = 1/pow(rhoa,16.0);
    t10 = 1/pow(3.0,0.333333333333333);
    t11 = 1/pow(3.141592653589793,6.666666666666667);
    t12 = pow(grada,10.0);
    t13 = 1/pow(rhoa,13.33333333333333);
    t14 = 1/pow(3.141592653589793,5.333333333333333);
    t15 = pow(grada,8.0);
    t16 = 1/pow(rhoa,10.66666666666667);
    t17 = 1/pow(3.141592653589793,4.0);
    t18 = pow(grada,6.0);
    t19 = 1/pow(rhoa,8.0);
    t20 = 1/pow(3.141592653589793,2.666666666666667);
    t21 = pow(grada,4.0);
    t22 = 1/pow(rhoa,5.333333333333333);
    t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
        t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
        t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
        t10*t11*t12*t13+1.0;
    t24 = pow(t23,0.024974);
    t25 = 1/pow(rhoa,2.333333333333333);
    t26 = 1/t5;
    t27 = pow(rhoa,0.333333333333333);
    t28 = 1/pow(rhoa,17.0);
    t29 = 1/pow(rhoa,14.33333333333333);
    t30 = 1/pow(rhoa,11.66666666666667);
    t31 = 1/pow(rhoa,9.0);
    t32 = 1/pow(rhoa,6.333333333333333);
    t33 = 1/pow(rhoa,3.666666666666667);
    t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
        t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
        t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
        t7*t8*t28;
    t35 = 1/pow(t23,0.975026);
    t36 = pow(rhoa,1.333333333333333);
    t37 = 1/t36;
    t38 = pow(grada,11.0);
    t39 = pow(grada,9.0);
    t40 = pow(grada,7.0);
    t41 = pow(grada,5.0);
    t42 = pow(grada,3.0);
    t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
        t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
        t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
        t7*t38*t9;
    t44 = 1/pow(t5,3.0);
    t45 = 1/pow(rhoa,6.0);
    t46 = 1/pow(rhoa,3.333333333333333);
    t47 = 1/pow(rhoa,0.666666666666667);
    t48 = pow(t34,2.0);
    t49 = 1/pow(t23,1.975026);
    t50 = 1/pow(rhoa,18.0);
    t51 = 1/pow(rhoa,15.33333333333333);
    t52 = 1/pow(rhoa,12.66666666666667);
    t53 = 1/pow(rhoa,10.0);
    t54 = 1/pow(rhoa,7.333333333333333);
    t55 = 1/pow(rhoa,4.666666666666667);
    t56 = 12.07834022939412*t1*t2*t3*t55+20.96333333333333*
        t10*t20*t21*t54+2.802125*t17*t18*t53+0.65457561728395*t1*t14*
        t15*t52+10.85299318415638*t10*t11*t12*t51+0.045864583333333*
        t7*t8*t50;
    t57 = 1/pow(rhoa,5.0);
    t58 = -6.588185579669522*t1*t2*grada*t33-13.24*t10*t20*
        t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*t1*
        t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*t7*
        t38*t28;
    t59 = 1/pow(rhoa,4.0);
    t60 = pow(t43,2.0);
    t61 = 2.470569592376071*t1*t2*t4+7.4475*t10*t20*t3*t22+
        1.167552083333333*t17*t21*t19+0.294559027777778*t1*t14*t18*
        t16+5.111002604166666*t10*t11*t15*t13+0.0222578125*t7*t12*
        t9;
    t62 = 1/pow(t5,4.0);
    t63 = 1/pow(t23,2.975026);

   /* code */
    res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
        t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
    res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
        t6*t24*t37;
    res[2] = 0.444444444444444*t24*t26*t47-6.666666666666668E-9*
        t1*t2*t3*t6*t24*t46+2.96296296296296E-17*t10*t20*t21*t44*t24*
        t45-0.024350299324*t48*t26*t49*t36+0.024974*t56*t26*t35*t36+
        0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*t1*
        t2*t3*t34*t6*t35*t25;
    res[3] = 0.024974*t58*t26*t35*t36-0.024350299324*t34*
        t43*t26*t49*t36+0.033298666666667*t43*t26*t35*t27-1.2487E-10*
        t1*t2*grada*t34*t6*t35*t37+6.666666666666667E-9*t1*t2*grada*
        t6*t24*t25+1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25-2.222222222222222E-17*
        t10*t20*t42*t44*t24*t57;
    res[4] = 0.024974*t61*t26*t35*t36-0.024350299324*t60*
        t26*t49*t36-5.0E-9*t1*t2*t6*t24*t37-2.4974E-10*t1*t2*grada*
        t43*t6*t35*t37+1.666666666666667E-17*t10*t20*t3*t44*t24*t59;
    res[5] = 0.024974*t26*t35*t36*(-0.8255625*t7*t8/pow(rhoa,
        19.0)-166.4125621570645*t10*t11*t12/pow(rhoa,16.33333333333333)-
        8.291291152263373*t1*t14*t15/pow(rhoa,13.66666666666667)-28.02125*
        t17*t18/pow(rhoa,11.0)-153.7311111111111*t10*t20*t21/pow(rhoa,
        8.333333333333334)-56.36558773717258*t1*t2*t3/pow(rhoa,5.666666666666667))+
        1.975308641975309E-25*t17*t18*t24*t62/pow(rhoa,9.666666666666666)-
        2.074074074074074E-16*t10*t20*t21*t24*t44/pow(rhoa,7.0)+2.518518518518519E-8*
        t1*t2*t24*t3*t6/pow(rhoa,4.333333333333333)-0.296296296296296*
        t24*t26/pow(rhoa,1.666666666666667)+0.048092474272682*t26*
        pow(t34,3.0)*t36*t63+0.033298666666667*t34*t26*t35*t47-4.9948E-10*
        t1*t2*t3*t34*t6*t35*t46+2.21991111111111E-18*t10*t20*t21*t34*
        t44*t35*t45-0.073050897972*t56*t34*t26*t49*t36-0.097401197296*
        t48*t26*t49*t27+0.099896*t56*t26*t35*t27-4.8700598648E-10*
        t1*t2*t3*t48*t6*t49*t25+4.9948E-10*t1*t2*t3*t56*t6*t35*t25;
    res[6] = -1.481481481481482E-25*t17*t24*t41*t62/pow(rhoa,
        8.666666666666666)-1.109955555555555E-18*t10*t20*t42*t34*t44*
        t35*t57+0.011099555555556*t43*t26*t35*t47-1.664933333333334E-10*
        t1*t2*t3*t43*t6*t35*t46-1.555555555555556E-8*t1*t2*grada*t6*
        t24*t46+7.3997037037037E-19*t10*t20*t21*t43*t44*t35*t45+1.407407407407407E-16*
        t10*t20*t42*t44*t24*t45+1.2175149662E-10*t1*t2*grada*t48*t6*
        t49*t37-1.2487E-10*t1*t2*grada*t56*t6*t35*t37+0.048092474272682*
        t48*t43*t26*t63*t36-0.024350299324*t56*t43*t26*t49*t36-0.048700598648*
        t58*t34*t26*t49*t36+0.024974*(24.15668045878825*t1*t2*grada*
        t55+83.85333333333332*t10*t20*t42*t54+16.81275*t17*t41*t53+
        5.236604938271604*t1*t14*t40*t52+108.5299318415638*t10*t11*
        t39*t51+0.550375*t7*t38*t50)*t26*t35*t36-0.064934131530667*
        t34*t43*t26*t49*t27+0.066597333333333*t58*t26*t35*t27-3.24670657653333E-10*
        t1*t2*t3*t34*t43*t6*t49*t25+3.32986666666667E-10*t1*t2*t3*
        t58*t6*t35*t25+3.32986666666667E-10*t1*t2*grada*t34*t6*t35*
        t25;
    res[7] = 1.111111111111111E-25*t17*t21*t24*t62/pow(rhoa,
        7.666666666666667)+4.16233333333333E-19*t10*t20*t3*t34*t44*
        t35*t59-1.109955555555555E-18*t10*t20*t42*t43*t44*t35*t57-
        8.88888888888889E-17*t10*t20*t3*t44*t24*t57+2.4350299324E-10*
        t1*t2*grada*t34*t43*t6*t49*t37-2.4974E-10*t1*t2*grada*t58*
        t6*t35*t37-1.2487E-10*t1*t2*t34*t6*t35*t37+0.048092474272682*
        t34*t60*t26*t63*t36-0.024350299324*t34*t61*t26*t49*t36-0.048700598648*
        t58*t43*t26*t49*t36+0.024974*(-6.588185579669522*t1*t2*t33-
        39.72*t10*t20*t3*t32-9.340416666666666*t17*t21*t31-3.141962962962963*
        t1*t14*t18*t30-68.14670138888889*t10*t11*t15*t29-0.356125*
        t7*t12*t28)*t26*t35*t36-0.032467065765333*t60*t26*t49*t27+
        0.033298666666667*t61*t26*t35*t27-1.623353288266666E-10*t1*
        t2*t3*t60*t6*t49*t25+1.664933333333333E-10*t1*t2*t3*t61*t6*
        t35*t25+3.32986666666667E-10*t1*t2*grada*t43*t6*t35*t25+6.666666666666667E-9*
        t1*t2*t6*t24*t25;
    res[8] = -8.33333333333333E-26*t17*t24*t42*t62/pow(rhoa,
        6.666666666666667)+0.048092474272682*t26*t36*pow(t43,3.0)*
        t63+1.2487E-18*t10*t20*t3*t43*t44*t35*t59+5.0E-17*t10*t20*
        grada*t44*t24*t59+3.6525448986E-10*t1*t2*grada*t60*t6*t49*
        t37-3.7461E-10*t1*t2*grada*t61*t6*t35*t37-3.7461E-10*t1*t2*
        t43*t6*t35*t37-0.073050897972*t61*t43*t26*t49*t36+0.024974*
        (14.895*t10*t20*grada*t22+4.670208333333333*t17*t42*t19+1.767354166666667*
        t1*t14*t41*t16+40.88802083333333*t10*t11*t40*t13+0.222578125*
        t7*t39*t9)*t26*t35*t36;

}

static void
lg93_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    lg93_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        lg93_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
lg93_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71, t72, t73, t74;
    real t75, t76, t77, t78, t79, t80, t81, t82;
    real t83, t84, t85;

    t1 = 1/pow(3.0,0.666666666666667);
    t2 = 1/pow(3.141592653589793,1.333333333333333);
    t3 = pow(grada,2.0);
    t4 = 1/pow(rhoa,2.666666666666667);
    t5 = 2.5E-9*t1*t2*t3*t4+1.0;
    t6 = 1/pow(t5,2.0);
    t7 = 1/pow(3.141592653589793,8.0);
    t8 = pow(grada,12.0);
    t9 = 1/pow(rhoa,16.0);
    t10 = 1/pow(3.0,0.333333333333333);
    t11 = 1/pow(3.141592653589793,6.666666666666667);
    t12 = pow(grada,10.0);
    t13 = 1/pow(rhoa,13.33333333333333);
    t14 = 1/pow(3.141592653589793,5.333333333333333);
    t15 = pow(grada,8.0);
    t16 = 1/pow(rhoa,10.66666666666667);
    t17 = 1/pow(3.141592653589793,4.0);
    t18 = pow(grada,6.0);
    t19 = 1/pow(rhoa,8.0);
    t20 = 1/pow(3.141592653589793,2.666666666666667);
    t21 = pow(grada,4.0);
    t22 = 1/pow(rhoa,5.333333333333333);
    t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
        t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
        t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
        t10*t11*t12*t13+1.0;
    t24 = pow(t23,0.024974);
    t25 = 1/pow(rhoa,2.333333333333333);
    t26 = 1/t5;
    t27 = pow(rhoa,0.333333333333333);
    t28 = 1/pow(rhoa,17.0);
    t29 = 1/pow(rhoa,14.33333333333333);
    t30 = 1/pow(rhoa,11.66666666666667);
    t31 = 1/pow(rhoa,9.0);
    t32 = 1/pow(rhoa,6.333333333333333);
    t33 = 1/pow(rhoa,3.666666666666667);
    t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
        t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
        t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
        t7*t8*t28;
    t35 = 1/pow(t23,0.975026);
    t36 = pow(rhoa,1.333333333333333);
    t37 = 1/t36;
    t38 = pow(grada,11.0);
    t39 = pow(grada,9.0);
    t40 = pow(grada,7.0);
    t41 = pow(grada,5.0);
    t42 = pow(grada,3.0);
    t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
        t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
        t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
        t7*t38*t9;
    t44 = 1/pow(t5,3.0);
    t45 = 1/pow(rhoa,6.0);
    t46 = 1/pow(rhoa,3.333333333333333);
    t47 = 1/pow(rhoa,0.666666666666667);
    t48 = pow(t34,2.0);
    t49 = 1/pow(t23,1.975026);
    t50 = 1/pow(rhoa,18.0);
    t51 = 1/pow(rhoa,15.33333333333333);
    t52 = 1/pow(rhoa,12.66666666666667);
    t53 = 1/pow(rhoa,10.0);
    t54 = 1/pow(rhoa,7.333333333333333);
    t55 = 1/pow(rhoa,4.666666666666667);
    t56 = 12.07834022939412*t1*t2*t3*t55+20.96333333333333*
        t10*t20*t21*t54+2.802125*t17*t18*t53+0.65457561728395*t1*t14*
        t15*t52+10.85299318415638*t10*t11*t12*t51+0.045864583333333*
        t7*t8*t50;
    t57 = 1/pow(rhoa,5.0);
    t58 = -6.588185579669522*t1*t2*grada*t33-13.24*t10*t20*
        t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*t1*
        t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*t7*
        t38*t28;
    t59 = 1/pow(rhoa,4.0);
    t60 = pow(t43,2.0);
    t61 = 2.470569592376071*t1*t2*t4+7.4475*t10*t20*t3*t22+
        1.167552083333333*t17*t21*t19+0.294559027777778*t1*t14*t18*
        t16+5.111002604166666*t10*t11*t15*t13+0.0222578125*t7*t12*
        t9;
    t62 = 1/pow(t5,4.0);
    t63 = 1/pow(rhoa,9.666666666666666);
    t64 = 1/pow(rhoa,7.0);
    t65 = 1/pow(rhoa,4.333333333333333);
    t66 = 1/pow(rhoa,1.666666666666667);
    t67 = pow(t34,3.0);
    t68 = 1/pow(t23,2.975026);
    t69 = 1/pow(rhoa,19.0);
    t70 = 1/pow(rhoa,16.33333333333333);
    t71 = 1/pow(rhoa,13.66666666666667);
    t72 = 1/pow(rhoa,11.0);
    t73 = 1/pow(rhoa,8.333333333333334);
    t74 = 1/pow(rhoa,5.666666666666667);
    t75 = -56.36558773717258*t1*t2*t3*t74-153.7311111111111*
        t10*t20*t21*t73-28.02125*t17*t18*t72-8.291291152263373*t1*
        t14*t15*t71-166.4125621570645*t10*t11*t12*t70-0.8255625*t7*
        t8*t69;
    t76 = 1/pow(rhoa,8.666666666666666);
    t77 = 24.15668045878825*t1*t2*grada*t55+83.85333333333332*
        t10*t20*t42*t54+16.81275*t17*t41*t53+5.236604938271604*t1*
        t14*t40*t52+108.5299318415638*t10*t11*t39*t51+0.550375*t7*
        t38*t50;
    t78 = 1/pow(rhoa,7.666666666666667);
    t79 = -6.588185579669522*t1*t2*t33-39.72*t10*t20*t3*t32-
        9.340416666666666*t17*t21*t31-3.141962962962963*t1*t14*t18*
        t30-68.14670138888889*t10*t11*t15*t29-0.356125*t7*t12*t28;
    t80 = 1/
        pow(rhoa,6.666666666666667);
    t81 = pow(t43,3.0);
    t82 = 14.895*t10*t20*grada*t22+4.670208333333333*t17*
        t42*t19+1.767354166666667*t1*t14*t41*t16+40.88802083333333*
        t10*t11*t40*t13+0.222578125*t7*t39*t9;
    t83 = 1/pow(t5,5.0);
    t84 = 1/pow(t23,3.975026);
    t85 = 1/pow(rhoa,9.333333333333334);

   /* code */
    res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
        t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
    res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
        t6*t24*t37;
    res[2] = 0.444444444444444*t24*t26*t47-6.666666666666668E-9*
        t1*t2*t3*t6*t24*t46+2.96296296296296E-17*t10*t20*t21*t44*t24*
        t45-0.024350299324*t48*t26*t49*t36+0.024974*t56*t26*t35*t36+
        0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*t1*
        t2*t3*t34*t6*t35*t25;
    res[3] = 0.024974*t58*t26*t35*t36-0.024350299324*t34*
        t43*t26*t49*t36+0.033298666666667*t43*t26*t35*t27-1.2487E-10*
        t1*t2*grada*t34*t6*t35*t37+6.666666666666667E-9*t1*t2*grada*
        t6*t24*t25+1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25-2.222222222222222E-17*
        t10*t20*t42*t44*t24*t57;
    res[4] = 0.024974*t61*t26*t35*t36-0.024350299324*t60*
        t26*t49*t36-5.0E-9*t1*t2*t6*t24*t37-2.4974E-10*t1*t2*grada*
        t43*t6*t35*t37+1.666666666666667E-17*t10*t20*t3*t44*t24*t59;
    res[5] = -0.296296296296296*t24*t26*t66+2.518518518518519E-8*
        t1*t2*t3*t6*t24*t65-2.074074074074074E-16*t10*t20*t21*t44*
        t24*t64+1.975308641975309E-25*t17*t18*t62*t24*t63+0.033298666666667*
        t34*t26*t35*t47-4.9948E-10*t1*t2*t3*t34*t6*t35*t46+2.21991111111111E-18*
        t10*t20*t21*t34*t44*t35*t45+0.048092474272682*t67*t26*t68*
        t36-0.073050897972*t56*t34*t26*t49*t36+0.024974*t75*t26*t35*
        t36-0.097401197296*t48*t26*t49*t27+0.099896*t56*t26*t35*t27-
        4.8700598648E-10*t1*t2*t3*t48*t6*t49*t25+4.9948E-10*t1*t2*
        t3*t56*t6*t35*t25;
    res[6] = 0.024974*t77*t26*t35*t36-0.024350299324*t56*
        t43*t26*t49*t36-0.048700598648*t58*t34*t26*t49*t36+0.048092474272682*
        t48*t43*t26*t68*t36+0.066597333333333*t58*t26*t35*t27-0.064934131530667*
        t34*t43*t26*t49*t27+0.011099555555556*t43*t26*t35*t47-1.2487E-10*
        t1*t2*grada*t56*t6*t35*t37+1.2175149662E-10*t1*t2*grada*t48*
        t6*t49*t37+3.32986666666667E-10*t1*t2*grada*t34*t6*t35*t25+
        3.32986666666667E-10*t1*t2*t3*t58*t6*t35*t25-3.24670657653333E-10*
        t1*t2*t3*t34*t43*t6*t49*t25-1.555555555555556E-8*t1*t2*grada*
        t6*t24*t46-1.664933333333334E-10*t1*t2*t3*t43*t6*t35*t46-1.109955555555555E-18*
        t10*t20*t42*t34*t44*t35*t57+1.407407407407407E-16*t10*t20*
        t42*t44*t24*t45+7.3997037037037E-19*t10*t20*t21*t43*t44*t35*
        t45-1.481481481481482E-25*t17*t41*t62*t24*t76;
    res[7] = 0.024974*t79*t26*t35*t36-0.048700598648*t58*
        t43*t26*t49*t36-0.024350299324*t34*t61*t26*t49*t36+0.048092474272682*
        t34*t60*t26*t68*t36+0.033298666666667*t61*t26*t35*t27-0.032467065765333*
        t60*t26*t49*t27-1.2487E-10*t1*t2*t34*t6*t35*t37-2.4974E-10*
        t1*t2*grada*t58*t6*t35*t37+2.4350299324E-10*t1*t2*grada*t34*
        t43*t6*t49*t37+6.666666666666667E-9*t1*t2*t6*t24*t25+3.32986666666667E-10*
        t1*t2*grada*t43*t6*t35*t25+1.664933333333333E-10*t1*t2*t3*
        t61*t6*t35*t25-1.623353288266666E-10*t1*t2*t3*t60*t6*t49*t25+
        4.16233333333333E-19*t10*t20*t3*t34*t44*t35*t59-8.88888888888889E-17*
        t10*t20*t3*t44*t24*t57-1.109955555555555E-18*t10*t20*t42*t43*
        t44*t35*t57+1.111111111111111E-25*t17*t21*t62*t24*t78;
    res[8] = 0.024974*t82*t26*t35*t36-0.073050897972*t61*
        t43*t26*t49*t36+0.048092474272682*t81*t26*t68*t36-3.7461E-10*
        t1*t2*t43*t6*t35*t37-3.7461E-10*t1*t2*grada*t61*t6*t35*t37+
        3.6525448986E-10*t1*t2*grada*t60*t6*t49*t37+5.0E-17*t10*t20*
        grada*t44*t24*t59+1.2487E-18*t10*t20*t3*t43*t44*t35*t59-8.33333333333333E-26*
        t17*t42*t62*t24*t80;
    res[9] = 0.024974*t26*t35*t36*(15.6856875*t7*t8/pow(rhoa,
        20.0)+2718.071848565387*t10*t11*t12/pow(rhoa,17.33333333333333)+
        113.3143124142661*t1*t14*t15/pow(rhoa,14.66666666666667)+308.23375*
        t17*t18/pow(rhoa,12.0)+1281.092592592593*t10*t20*t21*t85+319.4049971773113*
        t1*t2*t3*t80)-0.143076361365561*t26*pow(t34,4.0)*t36*t84-0.029598814814815*
        t34*t26*t35*t66+2.51589925925926E-9*t1*t2*t3*t34*t6*t35*t65-
        2.071917037037037E-17*t10*t20*t21*t34*t44*t35*t64+1.973254320987654E-26*
        t17*t18*t34*t62*t35*t63-0.073050897972*t26*t36*t49*pow(t56,
        2.0)-0.064934131530667*t48*t26*t49*t47+0.066597333333333*t56*
        t26*t35*t47+9.7401197296E-10*t1*t2*t3*t48*t6*t49*t46-9.9896E-10*
        t1*t2*t3*t56*t6*t35*t46-4.32894210204444E-18*t10*t20*t21*t48*
        t44*t49*t45+4.43982222222222E-18*t10*t20*t21*t56*t44*t35*t45+
        0.493827160493827*t24*t26*t4+0.288554845636095*t56*t48*t26*
        t68*t36-0.097401197296*t75*t34*t26*t49*t36+0.256493196120973*
        t67*t26*t68*t27-0.389604789184*t56*t34*t26*t49*t27+0.133194666666667*
        t75*t26*t35*t27+1.282465980604865E-9*t1*t2*t3*t67*t6*t68*t25-
        1.94802394592E-9*t1*t2*t3*t56*t34*t6*t49*t25+6.65973333333333E-10*
        t1*t2*t3*t75*t6*t35*t25-1.1111111111111112E-7*t1*t2*t3*t6*
        t24*t22+1.563786008230453E-15*t10*t20*t21*t44*t24*t19-3.29218106995885E-24*
        t17*t18*t62*t24*t16+5.26748971193416E-33*t1*t14*t15*t83*t24*
        t13;
    res[10] = -3.95061728395062E-33*t1*t14*t24*t40*t83/pow(rhoa,
        12.33333333333333)-1.109955555555556E-26*t17*t41*t34*t62*t35*
        t76-0.007399703703704*t43*t26*t35*t66+6.28974814814815E-10*
        t1*t2*t3*t43*t6*t35*t65+5.185185185185187E-8*t1*t2*grada*t6*
        t24*t65-5.17979259259259E-18*t10*t20*t21*t43*t44*t35*t64-9.1358024691358E-16*
        t10*t20*t42*t44*t24*t64+4.93313580246914E-27*t17*t18*t43*t62*
        t35*t63+2.22222222222222E-24*t17*t41*t62*t24*t63+1.623353288266667E-18*
        t10*t20*t42*t48*t44*t49*t57-1.664933333333333E-18*t10*t20*
        t42*t56*t44*t35*t57-0.032467065765333*t34*t43*t26*t49*t47+
        0.033298666666667*t58*t26*t35*t47+4.8700598648E-10*t1*t2*t3*
        t34*t43*t6*t49*t46-4.9948E-10*t1*t2*t3*t58*t6*t35*t46-1.165453333333334E-9*
        t1*t2*grada*t34*t6*t35*t46-2.16447105102222E-18*t10*t20*t21*
        t34*t43*t44*t49*t45+2.21991111111111E-18*t10*t20*t21*t58*t44*
        t35*t45+1.054457777777778E-17*t10*t20*t42*t34*t44*t35*t45-
        2.40462371363412E-10*t1*t2*grada*t67*t6*t68*t37+3.6525448986E-10*
        t1*t2*grada*t56*t34*t6*t49*t37-1.2487E-10*t1*t2*grada*t75*
        t6*t35*t37-0.143076361365561*t67*t43*t26*t84*t36+0.144277422818047*
        t58*t48*t26*t68*t36+0.144277422818047*t56*t34*t43*t26*t68*
        t36-0.073050897972*t56*t58*t26*t49*t36-0.024350299324*t75*
        t43*t26*t49*t36-0.073050897972*t77*t34*t26*t49*t36+0.024974*
        (-112.7311754743452*t1*t2*grada*t74-614.9244444444444*t10*
        t20*t42*t73-168.1275*t17*t41*t72-66.33032921810698*t1*t14*
        t40*t71-1664.125621570645*t10*t11*t39*t70-9.906749999999999*
        t7*t38*t69)*t26*t35*t36+0.19236989709073*t48*t43*t26*t68*t27-
        0.097401197296*t56*t43*t26*t49*t27-0.194802394592*t58*t34*
        t26*t49*t27+0.099896*t77*t26*t35*t27+9.61849485453648E-10*
        t1*t2*t3*t48*t43*t6*t68*t25-4.8700598648E-10*t1*t2*grada*t48*
        t6*t49*t25-4.8700598648E-10*t1*t2*t3*t56*t43*t6*t49*t25-9.7401197296E-10*
        t1*t2*t3*t58*t34*t6*t49*t25+4.9948E-10*t1*t2*t3*t77*t6*t35*
        t25+4.9948E-10*t1*t2*grada*t56*t6*t35*t25;
    res[11] = 2.962962962962963E-33*t1*t14*t18*t24*t83/pow(rhoa,
        11.33333333333333)+5.54977777777778E-27*t17*t21*t34*t62*t35*
        t78-7.3997037037037E-27*t17*t41*t43*t62*t35*t76-1.444444444444445E-24*
        t17*t21*t62*t24*t76-4.05838322066667E-19*t10*t20*t3*t48*t44*
        t49*t59+4.16233333333333E-19*t10*t20*t3*t56*t44*t35*t59-0.048700598648*
        t26*t36*t49*pow(t58,2.0)+2.16447105102222E-18*t10*t20*t42*
        t34*t43*t44*t49*t57-2.21991111111111E-18*t10*t20*t42*t58*t44*
        t35*t57-4.43982222222222E-18*t10*t20*t3*t34*t44*t35*t57-0.010822355255111*
        t60*t26*t49*t47+0.011099555555556*t61*t26*t35*t47+1.623353288266667E-10*
        t1*t2*t3*t60*t6*t49*t46-1.664933333333334E-10*t1*t2*t3*t61*
        t6*t35*t46-7.76968888888889E-10*t1*t2*grada*t43*t6*t35*t46-
        1.555555555555556E-8*t1*t2*t6*t24*t46-7.21490350340741E-19*
        t10*t20*t21*t60*t44*t49*t45+7.3997037037037E-19*t10*t20*t21*
        t61*t44*t35*t45+7.02971851851852E-18*t10*t20*t42*t43*t44*t35*
        t45+4.74074074074074E-16*t10*t20*t3*t44*t24*t45-4.80924742726824E-10*
        t1*t2*grada*t48*t43*t6*t68*t37+1.2175149662E-10*t1*t2*t48*
        t6*t49*t37+2.4350299324E-10*t1*t2*grada*t56*t43*t6*t49*t37+
        4.8700598648E-10*t1*t2*grada*t58*t34*t6*t49*t37-2.4974E-10*
        t1*t2*grada*t77*t6*t35*t37-1.2487E-10*t1*t2*t56*t6*t35*t37-
        0.143076361365561*t48*t60*t26*t84*t36+0.048092474272682*t48*
        t61*t26*t68*t36+0.048092474272682*t56*t60*t26*t68*t36+0.19236989709073*
        t58*t34*t43*t26*t68*t36-0.024350299324*t56*t61*t26*t49*t36-
        0.048700598648*t77*t43*t26*t49*t36-0.048700598648*t79*t34*
        t26*t49*t36+0.024974*(24.15668045878825*t1*t2*t55+251.56*t10*
        t20*t3*t54+84.06375*t17*t21*t53+36.65623456790123*t1*t14*t18*
        t52+976.769386574074*t10*t11*t15*t51+6.054125*t7*t12*t50)*
        t26*t35*t36+0.128246598060486*t34*t60*t26*t68*t27-0.064934131530667*
        t34*t61*t26*t49*t27-0.129868263061333*t58*t43*t26*t49*t27+
        0.066597333333333*t79*t26*t35*t27+6.41232990302432E-10*t1*
        t2*t3*t34*t60*t6*t68*t25-3.24670657653333E-10*t1*t2*t3*t34*
        t61*t6*t49*t25-6.49341315306667E-10*t1*t2*t3*t58*t43*t6*t49*
        t25-6.49341315306667E-10*t1*t2*grada*t34*t43*t6*t49*t25+3.32986666666667E-10*
        t1*t2*t3*t79*t6*t35*t25+6.65973333333333E-10*t1*t2*grada*t58*
        t6*t35*t25+3.32986666666667E-10*t1*t2*t34*t6*t35*t25;
    res[12] = -2.222222222222222E-33*t1*t14*t24*t41*t83/pow(rhoa,
        10.33333333333333)-2.08116666666667E-27*t17*t42*t34*t62*t35*
        t80+8.32466666666667E-27*t17*t21*t43*t62*t35*t78+8.88888888888889E-25*
        t17*t42*t62*t24*t78-1.2175149662E-18*t10*t20*t3*t34*t43*t44*
        t49*t59+1.2487E-18*t10*t20*t3*t58*t44*t35*t59+1.2487E-18*t10*
        t20*grada*t34*t44*t35*t59+1.623353288266667E-18*t10*t20*t42*
        t60*t44*t49*t57-1.664933333333333E-18*t10*t20*t42*t61*t44*
        t35*t57-6.65973333333333E-18*t10*t20*t3*t43*t44*t35*t57-2.0E-16*
        t10*t20*grada*t44*t24*t57-7.21387114090236E-10*t1*t2*grada*
        t34*t60*t6*t68*t37+3.6525448986E-10*t1*t2*grada*t34*t61*t6*
        t49*t37+7.3050897972E-10*t1*t2*grada*t58*t43*t6*t49*t37+3.6525448986E-10*
        t1*t2*t34*t43*t6*t49*t37-3.7461E-10*t1*t2*grada*t79*t6*t35*
        t37-3.7461E-10*t1*t2*t58*t6*t35*t37-0.143076361365561*t34*
        t81*t26*t84*t36+0.144277422818047*t58*t60*t26*t68*t36+0.144277422818047*
        t34*t61*t43*t26*t68*t36-0.073050897972*t58*t61*t26*t49*t36-
        0.073050897972*t79*t43*t26*t49*t36-0.024350299324*t82*t34*
        t26*t49*t36+0.024974*(-79.44*t10*t20*grada*t32-37.36166666666666*
        t17*t42*t31-18.85177777777778*t1*t14*t41*t30-545.1736111111111*
        t10*t11*t40*t29-3.56125*t7*t39*t28)*t26*t35*t36+0.064123299030243*
        t81*t26*t68*t27-0.097401197296*t61*t43*t26*t49*t27+0.033298666666667*
        t82*t26*t35*t27+3.20616495151216E-10*t1*t2*t3*t81*t6*t68*t25-
        4.8700598648E-10*t1*t2*grada*t60*t6*t49*t25-4.8700598648E-10*
        t1*t2*t3*t61*t43*t6*t49*t25+1.664933333333333E-10*t1*t2*t3*
        t82*t6*t35*t25+4.9948E-10*t1*t2*grada*t61*t6*t35*t25+4.9948E-10*
        t1*t2*t43*t6*t35*t25;
    res[13] = 1.666666666666667E-33*t1*t14*t21*t83*t24*t85-
        0.143076361365561*t26*t36*pow(t43,4.0)*t84-8.32466666666667E-27*
        t17*t42*t43*t62*t35*t80-5.0E-25*t17*t3*t62*t24*t80-0.073050897972*
        t26*t36*t49*pow(t61,2.0)-2.4350299324E-18*t10*t20*t3*t60*t44*
        t49*t59+2.4974E-18*t10*t20*t3*t61*t44*t35*t59+4.9948E-18*t10*
        t20*grada*t43*t44*t35*t59+5.0E-17*t10*t20*t44*t24*t59-9.61849485453648E-10*
        t1*t2*grada*t81*t6*t68*t37+7.3050897972E-10*t1*t2*t60*t6*t49*
        t37+1.46101795944E-9*t1*t2*grada*t61*t43*t6*t49*t37-4.9948E-10*
        t1*t2*grada*t82*t6*t35*t37-7.4922E-10*t1*t2*t61*t6*t35*t37+
        0.288554845636095*t61*t60*t26*t68*t36-0.097401197296*t82*t43*
        t26*t49*t36+0.024974*(14.895*t10*t20*t22+14.010625*t17*t3*
        t19+8.836770833333333*t1*t14*t21*t16+286.2161458333333*t10*
        t11*t18*t13+2.003203125*t7*t15*t9)*t26*t35*t36;

}

static void
lg93_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    lg93_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        lg93_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
